#####
## Replication File, CAD and AUS Analyses:
##
## "Toward a Theory and Measure of Racial Benevolence"
##
## Edana Beauvais and Adam M. Enders
##
## Published in Political Behavior, 2026
#####

# Loading libraries
library(car)
library(egg) 
library(ggplot2)
library(psych)
library(RColorBrewer)
library(tidyverse)

# Clearing environment
rm(list = ls())

# Set working directory!

# Loading dataset
data <- read.csv("racial_benev_CA.csv")

################################################################################

# RELIABILITY

myvars <- c("police", "adopt", "private_property", "jobs", "better_schools")

library(lattice)
# Scree plot
check.scree <- fa(data[myvars], fm = "pa", SMC = TRUE, rotate = "none")
check.scree$values/ncol(data[myvars])

tiff("scree-CA.tiff", width = 5, height = 5, units = "in", res=300, compression = 'lzw')
xyplot(check.scree$values ~ 1:ncol(data[myvars]),
       aspect = 1,
       type = "b",
       col = "black",
       xlab = "Factor",
       ylab = "Eigenvalue",
       pch = 16
)
dev.off()

# Exploratory factor analysis results
fac1 <- fa(cor(data[myvars]), nfactors = 1, fm = "pa", rotate = "none")
fac1

fac <- fa(cor(data[myvars]), nfactors = 2, fm = "pa", rotate = "none")
fac$loadings
fac$e.values

# Checking Cronbach's Alpha
psych::alpha(data[myvars]) 

# IRT for item and test information
library(mirt)
grm.mod <- mirt(data[myvars], 1, itemtype = "graded", SE = TRUE)
coef(grm.mod, IRTpars = TRUE)

M2(grm.mod, type = "C2")
itemfit(grm.mod)

plot(grm.mod, type = 'info')

tiff("info-CA.tiff", width = 6, height = 5, units = "in", res=300, compression = 'lzw')
plot(grm.mod, type = 'infotrace', main = "")
dev.off()

################################################################################

# WHO HAS BENEVOLENT RACIAL ATTITUDES?

# Proportion of respondents agreeing with each item
# Note that each item is coded so higher values = greater benevolence
# 5 = strongly agree, 4 = agree
prop.table(table(data$adopt)) 
prop.table(table(data$better_schools))
prop.table(table(data$private_property))
prop.table(table(data$jobs))
prop.table(table(data$police))

# Distribution of Indigenous resentment & Racial benevolence
hist1 <- histogram(~br_scale, 
                   data = data,
                   aspect = 1,
                   xlab = "Racial Benevolence",
                   ylim = c(-1, 30),
                   panel = function(...){
                     panel.histogram(...)
                     panel.abline(v=mean(data$br_scale), lty=2)
                   }
)

hist2 <- histogram(~ir_scale, 
                   data = data,
                   aspect = 1,
                   xlab = "Indigenous Resentment",
                   ylim = c(-1, 30),
                   panel = function(...){
                     panel.histogram(...)
                     panel.abline(v=mean(data$ir_scale), lty=2)
                   }
)

png("hist-CA.png", width = 7, height = 3.75, units = "in", res=300)
print(hist1, position = c(0, 0, .5, 1), more = TRUE)
print(hist2, position = c(.5, 0, 1, 1), more = FALSE)
dev.off()

################################################################################

# Standardizing variables:
data_std <- subset(data, select = c(w_id_scale, white_guilt, sd_scale, feel_indig, 
                                    spend_scale_r, spend_aidtopoor, spend_crime, spend_welfare,
                                    br_scale, ir_scale, college, age, income, female, quebec))

data_std <- as.data.frame(scale(data_std))
data_std$pid_recoded <- data$pid_recoded # Adding the non-std PID variable


# Standardized coefficients from a multivariate regression of both racial benevolence and Indigenous resentment on sociodemographic variables. Educational attainment, household income, gender, age, and region (South vs. non-South in the U.S., and Quebec vs. Canada outside of Quebec in Canada)
regmod <- lm(cbind(br_scale, ir_scale) ~ college + age + income + female + quebec, data = data_std)
summary(regmod)

v <- vcov(regmod)
co <- setNames(c(coef(regmod)), rownames(v))
linearHypothesis(NULL, c("br_scale:college = ir_scale:college"), coef. = co, vcov. = v) # p < 0.05
linearHypothesis(NULL, c("br_scale:age = ir_scale:age"), coef. = co, vcov. = v) # p < 0.10
linearHypothesis(NULL, c("br_scale:income = ir_scale:income"), coef. = co, vcov. = v) # Not sig 
linearHypothesis(NULL, c("br_scale:female = ir_scale:female"), coef. = co, vcov. = v) # p < 0.05
linearHypothesis(NULL, c("br_scale:quebec = ir_scale:quebec"), coef. = co, vcov. = v) # p <0.05

################################################################################

# VALIDITY
# Correlates: Resentment, feelings, social distance, White guilt, White ID
# Resentment
cor.test(data$br_scale, data$ir_scale)

# Feelings
cor.test(data$br_scale, data$feel_indig)
cor.test(data$ir_scale, data$feel_indig)

# Social distance
cor.test(data$br_scale, data$sd_scale)
cor.test(data$ir_scale, data$sd_scale)

# Spend
cor.test(data$br_scale, data$spend_scale)
cor.test(data$ir_scale, data$spend_scale)
# Reliability of spend:
spend_vars <- c("spend_welfare", "spend_env", "spend_science", "spend_foreignaid", "spend_schools", "spend_crime", "spend_childcare", "spend_aidtopoor")
psych::alpha(data[spend_vars], check.keys = T)

# White guilt
cor.test(data$br_scale, data$w_g_scale)
cor.test(data$ir_scale, data$w_g_scale)
# Reliability of white guilt:
guilt_vars <- c("guilt_1", "guilt_2", "guilt_3")
psych::alpha(data[guilt_vars]) 

# White ID
cor.test(data$br_scale, data$w_id_scale)
cor.test(data$ir_scale, data$w_id_scale)
# reliability White ID
whiteid_vars <- c("wid_whiteid", "wid_pride", "wid_common") # 0.72
psych::alpha(data[whiteid_vars]) 

# Regressions of resentment & benevolent racial attitudes on partisanship, age, educational attainment, household income, gender, and region. Standardized coefficients
mod1 <- lm(feel_indig ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
mod2 <- lm(sd_scale ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
mod3 <- lm(white_guilt ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
mod4 <- lm(w_id_scale ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
summary(mod4)

library(stargazer)
stargazer(mod1, mod2, mod3, mod4,
          align=TRUE, 
          covariate.labels=c("Racial Benevolence","Indigenous resentment", 
                             "PID (Cons)", "PID (NDP)", "PID (Other)", "PID (None/DK)",
                             "College", "Age", "Income", "Gender", "Quebec"), 
          no.space=TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          df = FALSE)

################################################################################

# POTENTIAL POLICY CONSEQUENCES
# Spending scale needs to be reverse coded so that it's in the same direction as US data
data$spend_scale_r <- max(data$spend_scale) + min(data$spend_scale) - data$spend_scale
data$spend_aidtopoor <- max(data$spend_aidtopoor) + min(data$spend_aidtopoor) - data$spend_aidtopoor
data$spend_welfare <- max(data$spend_welfare) + min(data$spend_welfare) - data$spend_welfare
data$spend_crime <- max(data$spend_crime) + min(data$spend_crime) - data$spend_crime


# Standardized coefficients from OLS regressions of government spending preferences on racial benevolence, Indigenous resentment, and controls.
m_spend1 <- lm(spend_scale_r ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
m_spend2 <- lm(spend_aidtopoor ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
m_spend3 <- lm(spend_welfare ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)
m_spend4 <- lm(spend_crime ~ br_scale + ir_scale + pid_recoded + college + age + income + female + quebec, data = data_std)

stargazer(m_spend1, m_spend2, m_spend3, m_spend4,
          align=TRUE, 
          covariate.labels=c("Racial Benevolence","Indigenous resentment", 
                             "PID (Cons)", "PID (NDP)", "PID (Other)", "PID (None/DK)",
                             "College", "Age", "Income", "Gender", "Quebec"), 
          no.space=TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          df = FALSE)

# Supplementary spending analysis: second additive model with controls for White guilt and White identity and a third interactive model with interactions between 1) racial benevolence and White guilt and 2) racial/Indigenous resentment and White identity.
summary(lm(spend_scale_r ~ br_scale + ir_scale + pid_recoded + ideo + college + age + income + female + quebec + w_id_scale, data = data)) 
summary(lm(spend_scale_r ~ br_scale + ir_scale + pid_recoded + ideo + college + age + income + female + quebec + white_guilt, data = data))
summary(lm(spend_scale_r ~ br_scale + ir_scale + pid_recoded + ideo + college + age + income + female + quebec + w_id_scale + white_guilt + br_scale*white_guilt + ir_scale*w_id_scale, data = data))

################################################################################

# Donation item 
library(nnet)
table(data$donation)
# Change the reference category to "Male"
# Changing reference category
data <- data %>% mutate(donation = factor(donation, levels = c("Empower", "Assimilate", "None")))

mod_donate <- multinom(donation ~ br_scale + ir_scale + pid_recoded + ideo + college + age + income + female + region_factor, Hess = TRUE, data = data)

stargazer(mod_donate)

library(ggeffects)
library(dplyr)
library(effects)

tiff("donateBR-CA.tiff", width = 7, height = 3.5, units = "in", res=300, compression = 'lzw')
mod_donate %>% ggeffect(terms = "br_scale")%>% # use [all] to get a smooth result
  ggplot(., aes(x, predicted), facet = TRUE) + 
  geom_line(aes(color = response.level)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = response.level), alpha = .25) +
  # facet_grid(.~predicted) +
  facet_wrap(. ~ response.level, dir = "h") +
  ylab("Predicted Probability") +
  xlab("Racial Benevolence") +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5),
        panel.grid.minor.x=element_blank(),
        legend.position = "none",
        aspect.ratio = 1,
        axis.line = element_line(size=1, colour = "black"),
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size=12),
        strip.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 12, face="bold"),
        axis.title.x = element_text(size = 12, face="bold"))
dev.off()

tiff("donateRR-CA.tiff", width = 7, height = 3.5, units = "in", res=300, compression = 'lzw')
mod_donate %>% 
  ggeffect(terms = "ir_scale") %>% #use [all] to get a smooth result
  ggplot(., aes(x, predicted), facet = TRUE) + 
  geom_line(aes(color = response.level)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = response.level), alpha = .25) +
  # facet_grid(.~predicted) +
  facet_wrap(. ~ response.level, dir = "h") +
  ylab("Predicted Probability") +
  xlab("Indigenous Resentment") +
  # scale_x_discrete(labels=c("empower" = "Educational\nEmpowerment", "assimilate" = "Educational\nImmersion\n(Assimilative Policy)",
  # "none" = "None")) +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5),
        panel.grid.minor.x=element_blank(),
        legend.position = "none",
        aspect.ratio = 1,
        axis.line = element_line(size=1, colour = "black"),
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size=12),
        strip.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 12, face="bold"),
        axis.title.x = element_text(size = 12, face="bold"))
dev.off()

################################################################################

# Extension & Replication - Australia
library(readstata13)
d <- read.dta13("aus-data.dta")
dim(d)
head(d)
d<- na.omit(d)

# Aggregate measure of spend, higher values = more spend
dat_spend <- subset(d, select = c("spend_gw", "spend_healthcare", "spend_police", "spend_welfare", "spend_pension"))
d$spend_scale <- apply(dat_spend, MARGIN = 1, FUN = mean)
summary(d$spend_scale)

# Reliability (Australia) 
myvars <- c("police", "adopt", "property", "jobs", "better_schools")

library(lattice)
library(psych)
# Scree plot (Australia) 
check.scree <- fa(d[myvars], fm = "pa", SMC = TRUE, rotate = "none")
check.scree$values/ncol(d[myvars])

tiff("scree-AU.tiff", width = 5, height = 5, units = "in", res=300, compression = 'lzw')
xyplot(check.scree$values ~ 1:ncol(d[myvars]),
       aspect = 1,
       type = "b",
       col = "black",
       xlab = "Factor",
       ylab = "Eigenvalue",
       pch = 16
)

dev.off()

# Exploratory factor analysis results (Australia) 
fac1 <- fa(cor(d[myvars]), nfactors = 1, fm = "pa", rotate = "none")
fac1

fac <- fa(cor(d[myvars]), nfactors = 2, fm = "pa", rotate = "none")
fac$loadings
fac$e.values

# Checking Cronbach's Alpha (Australia) 
psych::alpha(d[myvars]) 

# IRT for item and test information (Australia) 
library(mirt)
grm.mod <- mirt(d[myvars], 1, itemtype = "graded", SE = TRUE)
coef(grm.mod, IRTpars = TRUE)

M2(grm.mod, type = "C2")
itemfit(grm.mod)

plot(grm.mod, type = 'info')

tiff("info-AU.tiff", width = 6, height = 5, units = "in", res=300, compression = 'lzw')
plot(grm.mod, type = 'infotrace', main = "")
dev.off()

irt.scores <- fscores(grm.mod, full.scores.SE = TRUE)
empirical_rxx(irt.scores)

# Distributions (Australia) 
hist1 <- histogram(~br_scale_std, 
                   data = dat,
                   aspect = 1,
                   xlab = "Racial Benevolence (Aus)",
                   ylim = c(-1, 30),
                   panel = function(...){
                     panel.histogram(...)
                     panel.abline(v=mean(data$br_scale), lty=2)
                   }
)

hist2 <- histogram(~ir_scale_std, 
                   data = dat,
                   aspect = 1,
                   xlab = "Indigenous Resentment (Aus)",
                   ylim = c(-1, 30),
                   panel = function(...){
                     panel.histogram(...)
                     panel.abline(v=mean(data$ir_scale), lty=2)
                   }
)

png("hist-AU.png", width = 7, height = 3.75, units = "in", res=300)
print(hist1, position = c(0, 0, .5, 1), more = TRUE)
print(hist2, position = c(.5, 0, 1, 1), more = FALSE)
dev.off()


# Validity (Australia) 
cor.test(d$br_scale, d$ir_scale)
# ideology (Australia) 
cor.test(d$br_scale, d$ideol)
cor.test(d$ir_scale, d$ideol)

# Spend (Australia) 
cor.test(d$br_scale, d$spend_scale) 
cor.test(d$ir_scale, d$spend_scale) 

# social distance (Australia) 
cor.test(d$br_scale, d$no_marry) 
cor.test(d$ir_scale, d$no_marry) 

# Guilt (Australia) 
cor.test(d$br_scale, d$guilt_scale) 
cor.test(d$ir_scale_std, d$guilt_scale) 

# White ID (Australia) 
cor.test(d$br_scale, d$white_id) 
cor.test(d$ir_scale, d$white_id) 

library(reshape2)
library(rms)

# Donation outcome (Australia) 

mod_donate <- multinom(donate ~ as.numeric(br_scale) + as.numeric(ir_scale) + as.factor(pid) + ideol + college + age_cat + female + income + region_dummy, Hess = TRUE, data = d)

stargazer(mod_donate)

library(ggeffects)
library(dplyr)
library(effects)
library(broom)
library(tidyr)
library(broom)

################################################################################

# CES comparison for appendix

dm <- read.csv("ces_compare.csv")
library(ggplot2)
library(ggthemes)
View(dm)
dm$var <- reorder(dm$var, dm$est)

tiff("compare-CA.tiff", width = 5, height = 5, units = "in", res=300, compression = 'lzw')
ggplot(dm, aes(y = est, x = var, ymin = lb, ymax = ub, color=as.factor(survey))) +
  geom_pointrange(size = .2) +
  coord_flip() +
  theme_bw() +
  scale_x_discrete("Variable") +
  scale_y_continuous("Mean", limits = c(0, 1)) +
  theme(aspect.ratio = 1.65, legend.title=element_blank())
dev.off()


